myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = mean.R1), , colour=I("gray50"), size=I(.06)) + 
scale_fill_gradientn(colors = myelin.colorbar, limits = c(0.54, 0.58), oob = squish, na.value = "gray75") + theme(legend.position = "none")

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/R1_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.5, height = 3)

myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = MWF), colour=I("gray50"), size=I(.06)) + 
scale_fill_gradientn(colors = myelin.colorbar, limits = c(0.02, 0.05), oob = squish, na.value = "gray75") + theme(legend.position = "none")

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/MWF_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.5, height = 3)

cor.test(myelin.maps$mean.R1, myelin.maps$MWF, method = c("spearman"))

myelin.maps %>% filter(MWF < 0.08) %>% #remove 2 MWF major outliers in PHA1
ggplot(., aes(x = MWF, y = mean.R1)) +
  geom_point(color = c("#8AABD5FF"), size = 0.2) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  labs(x="\nMyelin Water Fraction", y="R1\n") +
  scale_y_continuous(limits = c(0.525, 0.61), breaks = c(0.53, 0.55, 0.57, 0.59, 0.61)) +
  theme(
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/MWF_R1_correlationplot.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)

myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = MBP), colour=I("gray50"), size=I(.06)) + 
scale_fill_gradientn(colors = myelin.colorbar, limits = c(-0.65, 1), oob = squish, na.value = "gray75") + theme(legend.position = "none")

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/MBP_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.5, height = 3)

cor.test(myelin.maps$mean.R1, myelin.maps$MBP, method = c("spearman"))

myelin.maps %>% 
ggplot(., aes(x = MBP, y = mean.R1)) +
  geom_point(color = c("#8AABD5FF"), size = 0.2) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  labs(x="\nMyelin Basic Protein Expression (z)", y="R1\n") +
  scale_y_continuous(limits = c(0.525, 0.61), breaks = c(0.53, 0.55, 0.57, 0.59, 0.61)) +
  theme(
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/MBP_R1_correlationplot.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)

myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = SA.axis), colour=I("gray50"), size=I(.06)) + 
scale_fill_gradientn(colors = myelin.colorbar, trans = 'reverse', na.value = "gray75") + theme(legend.position = "none")

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/SAaxis_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.5, height = 3)

cor.test(myelin.maps$mean.R1, myelin.maps$SA.axis, method = c("spearman"))

myelin.maps %>%
ggplot(., aes(x = SA.axis, y = mean.R1)) +
  geom_point(color = c("#8AABD5FF"), size = 0.2) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  labs(x="\nSensorimotor-Association Axis", y="R1\n") +
  scale_y_continuous(limits = c(0.525, 0.61), breaks = c(0.53, 0.55, 0.57, 0.59, 0.61)) +
  theme(
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/SAaxis_R1_correlationplot.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)

#Compute average frontal cortex volume fraction at each depth for final study participants
cortexpve <- read_parquet("/Volumes/Hera/Projects/corticalmyelin_development/BIDS/derivatives/surface_metrics/7T_BrainMechanisms_cortexpve.parquet") #parquet with cortex fraction measures
cortexpve <- cortexpve %>% filter(atlas == "glasser") #just for glasser atlas 
cortexpve <- cortexpve[!(cortexpve$StructName %in% glasser.snr.exclude$orig_parcelname),] #just for good regions
cortexpve <- cortexpve[cortexpve$StructName %in% glasser.frontal$orig_parcelname,] #just for frontal lobe

extract_depth_cortexpve <- function(mymeasure){
  measure.df <- cortexpve %>% select(subject_id, session_id, StructName, all_of(mymeasure))
  names(measure.df) <- c("subject_id", "session_id", "orig_parcelname", "cortex.fraction")
  depth.cortexpve <- merge(measure.df, participants, by = c("subject_id", "session_id")) %>% select(cortex.fraction) %>% colMeans(.$cortex.fraction) #get cortex.fraction for final study participants only, and then average for this depth
  return(depth.cortexpve)}

pve.measures <- list("Mean_cortexpve.1.00%", "Mean_cortexpve.0.9%", "Mean_cortexpve.0.8%","Mean_cortexpve.0.7%","Mean_cortexpve.0.6%","Mean_cortexpve.0.5%","Mean_cortexpve.0.4%","Mean_cortexpve.0.3%","Mean_cortexpve.0.2%","Mean_cortexpve.0.1%", "Mean_cortexpve.0.0%") #measures to extract data for 
cortexpve.frontallobe <- lapply(pve.measures, function(x) { 
  extract_depth_cortexpve(x)}) 
cortexpve.frontallobe <- do.call(rbind, cortexpve.frontallobe) %>% as.data.frame()
cortexpve.frontallobe$depth <- c("a", "b", "c", "d", "e", "f", "g", "h", "i", "j", "k")
cortexpve.frontallobe$depth <- factor(cortexpve.frontallobe$depth, ordered = T, levels = c("k", "j", "i", "h", "g", "f", "e", "d", "c", "b", "a"))
kable(cortexpve.frontallobe)

cortexpve.frontallobe %>% 
  ggplot(aes(x = cortex.fraction, y = depth)) +
  geom_point(color = "#809dc4", shape = 15, size = 2) +
  theme_classic() +
  xlab("\nCortex Fraction") +
  ylab("Cortical Depth\n") +
  theme(
  legend.position = "none",
  axis.text.x = element_text(size = 6, family = "Arial", color = c("black")),
  axis.text.y = element_blank(),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks.x = element_line(linewidth = .2),
  axis.ticks.y = element_blank()) 

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1_supplementary/Depthplot_cortexfraction_frontal.pdf", device = "pdf", dpi = 500, width = 2.5, height = 2.5)

#Compute mean R1 in each region at each cortical depth
regional_means <- function(input.df){
  meanmap <- input.df %>% select(contains("ROI")) %>% colMeans() %>% as.data.frame() %>% set_names("mean.R1")
  meanmap$orig_parcelname <- row.names(meanmap)
  meanmap <- merge(meanmap, glasser.regions, by = c("orig_parcelname"), sort = F)
  return(meanmap)
}

regional.myelin.alldepths <- lapply(myelin.glasser.7T, function(depth) {
  regional_means(depth)
})

regional.myelin.alldepths <-  do.call(rbind, regional.myelin.alldepths) #regional mean R1 at all depths
regional.myelin.alldepths$depth <- substr(row.names(regional.myelin.alldepths), 1, 7) #assign depth variable
regional.myelin.alldepths$depth <- factor(regional.myelin.alldepths$depth, ordered = T, levels = ordered_depths)

regional.myelin.alldepths <-  regional.myelin.alldepths[!(regional.myelin.alldepths$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),] #exclude low SNR parcels
regional.myelin.alldepths <- regional.myelin.alldepths[regional.myelin.alldepths$orig_parcelname %in% glasser.frontal$orig_parcelname,] #only studying frontal


regional.myelin.alldepths %>% 
  ggplot(aes(x = mean.R1, y = depth)) +
  geom_violin(aes(fill = depth), color = "white") +
  stat_summary(fun = median, geom = "point", shape = 23, size = 2.25, color = "white") +
  theme_classic() +
  scale_fill_manual(values = depth_colorbar) +
  scale_x_continuous(limits = c(0.49, 0.65), expand = c(0,0)) +
  scale_y_discrete(expand = c(0, 0)) +
  xlab("\nR1") +
  ylab("Cortical Depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/Depthplot_regionalR1_frontal.pdf", device = "pdf", dpi = 500, width = 2.6, height = 4.9)

for(this.depth in c(unique(regional.myelin.alldepths$depth))){
  
  depth.data <- regional.myelin.alldepths %>% filter(depth == this.depth) %>% filter(label != "lh_L_10pp")
  
  depth.plot <- ggplot() + 
  geom_brain(data = depth.data, atlas = glasser,  mapping = aes(fill = mean.R1), colour=I("gray10"), size=I(.08)) + theme_classic() + theme(legend.position = "none") + 
   paletteer::scale_fill_paletteer_c("grDevices::PuBuGn", direction = -1, limits = c(0.48, 0.61), oob = squish, na.value = "white") + theme(legend.position = "none") +
  theme(legend.position = "none", axis.text = element_blank(), axis.line = element_blank(), axis.ticks = element_blank()) +
  new_scale_fill() + 
  geom_brain(data = glasser.plotting, atlas = glasser, mapping = aes(fill = cortex),  colour=I(alpha("gray10", 0)), size=I(0)) +
  scale_fill_manual(values = c(alpha("white", 0), "gray75")) 

  print(depth.plot)
  ggsave(filename = sprintf("/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/%s_corticalmap_frontal.pdf", this.depth), device = "pdf", dpi = 500, width = 4.95, height = 2.25)
}

examplar.areas.alldepth <- rbind(area4.myelin.alldepths, area46.myelin.alldepths, areaa24.myelin.alldepths)

ggplot() +
  geom_smooth(data = examplar.areas.alldepth, aes(x = mean.R1, y = depth, group = orig_parcelname, color= orig_parcelname), se = F, linewidth = .8) + 
  scale_color_manual(values = c(myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2], myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2])) +
  theme_classic() +
  xlab("\nR1") +
  ylab("Cortical Depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/Area_profiles_group.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)

#R1 values in all regions at all depths for all participants
myelin.glasser.7T.alldepths <-  do.call(rbind, myelin.glasser.7T)
myelin.glasser.7T.alldepths$depth <- substr(row.names(myelin.glasser.7T.alldepths), 1, 7) #assign depth variable
myelin.glasser.7T.alldepths$depth <- factor(myelin.glasser.7T.alldepths$depth, ordered = T, levels = ordered_depths)

#Wide to long format
cols_to_pivot <- names(myelin.glasser.7T.alldepths)[grepl("ROI", names(myelin.glasser.7T.alldepths))]
myelin.glasser.7T.alldepths <- myelin.glasser.7T.alldepths %>% pivot_longer(cols = all_of(cols_to_pivot), names_to = "orig_parcelname", values_to = "R1")

#Select exemplar regions of interest
myelin.glasser.7T.alldepths <- myelin.glasser.7T.alldepths %>% filter(orig_parcelname == "L_4_ROI" | orig_parcelname == "R_4_ROI" | orig_parcelname == "L_46_ROI" | orig_parcelname == "R_46_ROI" | orig_parcelname == "L_a24_ROI" | orig_parcelname == "R_a24_ROI")

#Format grouping variables
myelin.glasser.7T.alldepths$subject_id <- as.factor(myelin.glasser.7T.alldepths$subject_id)
myelin.glasser.7T.alldepths$orig_parcelname <- as.factor(myelin.glasser.7T.alldepths$orig_parcelname)

#Plot for the oldest 3 participants 
myelin.glasser.7T.alldepths <- myelin.glasser.7T.alldepths %>% filter(age > 32)

ggplot() +
  geom_smooth(data = myelin.glasser.7T.alldepths, aes(x = R1, y = depth, color = orig_parcelname, group = interaction(orig_parcelname, subject_id)), se = F, linewidth = .4) + 
  scale_color_manual(values = c(myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2], myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2])) +
  theme_classic() +
  xlab("\nR1") +
  ylab("Cortical Depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/Area_profiles_individual.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)

regional.myelin.alldepths.wide <- regional.myelin.alldepths %>% pivot_wider(id_cols = orig_parcelname, names_from = depth, values_from = mean.R1) %>% select(-orig_parcelname)
regional.myelin.alldepths.wide <- regional.myelin.alldepths.wide %>% select(depth_7, depth_6, depth_5, depth_4, depth_3, depth_2, depth_1)
depth.R1.cors <- cor(regional.myelin.alldepths.wide) # compute correlation matrix

ggcorrplot(corr = depth.R1.cors, method = "square", type = "upper", show.diag = T, lab = F, lab_size = ) +
scale_fill_gradientn(colors = c("#edf6ff", "#c2d4ed", "#97adcc", "#345c91"))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/Depth_correlationmatrix_frontal.pdf", device = "pdf", dpi = 500, width = 4.4, height = 2.2)